% This script generates ensembles of temperature trajectories for figure 2 in the SI.
close all
clear all
 
% Load RCP data 
load rf_rcp8_5.txt
yrs=rf_rcp8_5(:,1);
forcing=rf_rcp8_5(:,2);

% Set start time in 1850
tstart=86;
tend=252;

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Specify initial temperature anomaly in 1850
DeltaT0 =0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Specify model parameters %%%%%%%%%%%%%%%%%%
% We will measure time in years,
% Specify size of  time interval between points.
% so need a factor of number of seconds in year
% Effectively divides time in seconds by yearlength
dt=1; 
% So then have to multiply speed by yearlength
yearlength=3.1 *1e7;

% Specify number of points at which solution is wanted
nPeriods=length(forcing(tstart:end)); % Same as RCP dataset

% Specify # of intermediate steps taken between the sample points
nSteps=12; % Monthly

% Number of temperature trajectories without noise
Ntrials0=1;

% Number of temperature trajectories with noise
Ntrials=100; % Generating many small ensembles from which to estimate the time series properties of temperature trajectories, and match them to HadCRUT.

% Specify set of lambas to be simulated
ECS = 0.5:0.25:10;
lambdas = 3.7 ./ ECS;

% Specify set of noise variance to be simulated
sigmaQs = (0.2:0.1:2.0)*1e8;

% Choose set of effective heat capacities to be simulated
Ceffs = (0.2:0.2:2)*1e9;

parfor k = 1:1:length(lambdas)
    CS = ECS(k);
    lambda=lambdas(k);   % Joules/metres^2/Kelvin
    % And convert F/lambda to a function
    Level=ts2func(forcing(tstart:tend)/lambda);

    for j = 1:1:length(sigmaQs)
        sigmaQ = sigmaQs(j);
        for i = 1:1:length(Ceffs)
                Ceff=Ceffs(i)
                Sigma=sigmaQ/Ceff;
                Speed = lambda*yearlength/Ceff; % Kelvin per year 

                [DeltaT,Times]=Hasselmann(Speed,Level,Sigma,DeltaT0,Ntrials,nPeriods,nSteps,dt);
                %Create a CSV write file
                % First suppress the unit 3rd dimension in Delta T matrix
                Temp=squeeze(DeltaT);
                % Then concatenate 
                M=[Times Temp];
 
                filename = ['Historical_' num2str(CS) '_' num2str(sigmaQ) '_' num2str(Ceff) '_.txt'];
                csvwrite(filename,M)
        end
    end
end